\chapter{Implementation details}
Below, we provide three principal steps of the lattice Boltzmann simulation as examples of the C++  code implementation. The code shown in this chapter is from the 2D Single Phase LB simulator.\\


The collision step, described in Equation {\ref{collision_law}} can be implemented as a subroutine: $collision()$

\begin{lstlisting}

void collision() //particle collision incorporated into PDFs, if solid then bounceback boundary is used
{

// f[i][j][k] is the kth value of the distribution function of the particle on (x,y)

double GuoF,lm0,lm1;
for(int i=0;i<=NX;i++)	
	for(int j=0;j<=NY;j++)
	if (!Solid[i][j])// NOT INTERIOR SOLID NODE     
        	for(int k=0;k<Q;k++)
		{       

		//BODY FORCE CALCULATION, GuoF is a variable to incorporate a body force		lm0=((e[k][0]-u[i][j][0])*forcex[i][j]+(e[k][1]-u[i][j][1])*forcey[i][j])*3;
		lm1=(e[k][0]*u[i][j][0]+e[k][1]*u[i][j][1])*(e[k][0]*forcex[i][j]+e[k][1]*forcey[i][j])*9;
		GuoF=(1-1/(2*tau_f))*w[k]*(lm0+lm1);


		//The last term is an external body force term

		f[i][j][k]=f[i][j][k]+(feq(i,j,k,rho[i][j],u[i][j])-f[i][j][k])/tau_f+tau_f*GuoF;
		}
                else    
		{standard_bounceback_boundary(i,j);};    //BOUNCEBACK BOUNDARY CONDITION
		
}

\end{lstlisting}

The streaming step (Equation \ref{streaming_step_equ}) is implemented as a subroutine: $streaming()$
\begin{lstlisting}

void streaming()
{
int ip,jp;
for(int i=0;i<=NX;i++)	
for(int j=0;j<=NY;j++)
	if (!Solid[i][j])
        for(int k=0;k<Q;k++)
	{       
	ip=i-e[k][0];if (ip<0) {ip=NX;}; if (ip>NX) {ip=0;};
	jp=j-e[k][1];if (jp<0) {jp=NY;}; if (jp>NY) {jp=0;};
		
	F[i][j][k]=f[ip][jp][k];
		
	};

for(int i=0;i<=NX;i++)	
for(int j=0;j<=NY;j++)
	{                 
        f[i][j][0]=F[i][j][0];f[i][j][1]=F[i][j][1];f[i][j][2]=F[i][j][2];
	f[i][j][3]=F[i][j][3];f[i][j][4]=F[i][j][4];f[i][j][5]=F[i][j][5];
	f[i][j][6]=F[i][j][6];f[i][j][7]=F[i][j][7];f[i][j][8]=F[i][j][8];
	};
}
\end{lstlisting}

The macroscopic quantities including density $\rho$, velocity $v_x,v_y,v_z$ are calculated using  Equations (\ref{macro_cal_1},\ref{macro_cal_2}) in subroutine $comput\_macro\_variables()$:

\begin{lstlisting}
void comput_macro_variables()
{
for(int i=0;i<=NX;i++)	
for(int j=0;j<=NY;j++)
{ 

	if (!Solid[i][j]) // NOT INTERIOR SOLID NODE
	{
	u0[i][j][0]=u[i][j][0];
	u0[i][j][1]=u[i][j][1];
	rho[i][j]=0;
	u[i][j][0]=0;
	u[i][j][1]=0;
	for(int k=0;k<Q;k++)		{		
		rho[i][j]+=f[i][j][k];
		u[i][j][0]+=e[k][0]*f[i][j][k];
		u[i][j][1]+=e[k][1]*f[i][j][k];
		}
		
		
	u[i][j][0]=(u[i][j][0]+tau_f*forcex[i][j]/2)/rho[i][j];
	u[i][j][1]=(u[i][j][1]+tau_f*forcey[i][j]/2)/rho[i][j];

	}
                       
}

}

\end{lstlisting}
